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We study the influence of misalignment of the ferromagnetic exchange field on the equilibrium 
properties of hybrid structures, composed of superconducting (S) and ferromagnetic (F) parts. In 
particular, we study numerically the superconducting critical temperature T c in F-S-F trilayers and 
in F-S-F-S-F Josephson junctions as a function of the misalignment angle 9 of the ferromagnetic 
magnetization. We discuss the corresponding phase diagrams for these hybrid structures. For the 
Josephson junctions, a transition between the zero-phase and the 7r-phase ground state as a function 
of 9 takes place under certain conditions. Within the quasiclassical Green's function technique in 
the diffusive limit, we introduce a fast and effective method for calculating T c in such multilayer 
structures. 



I. INTRODUCTION 

The interest in superconductor-ferromagnet (S-F) hy- 
brid structures has considerably increased in the last 
decade due to their relevance for the development of 
nanometer scale electronic devices. The understanding 
of the superconducting proximity effect in S-F devices is 
of vital importance for such a goal. Consequently, ex- 
perimental and theoretical studies have focused on the 
influence that proximity induced spin-triplet pairing am- 
plitudes in S-F hybrid structures have on superconduct- 
ing properties of the entire structure. Among those are 
for example changes of the superconducting transition 
temperature T c of the device, or the switching between 
0-junctions and 7r-junctions as ground states in S-F-S 
Josephson devices as a function of some control parame- 
ter. 

The superconducting critical temperature T c 
in diffusive hybrid S-F structures has been 
studied both theoretically 1 ! 2 ! 3 ! 4 ! 5 ! 6 ! 7 ! 8 ! 9 ! 10 ! 11 and 
experimentally 12 ! 13 ! 14 ! 15 ! 16 ! 17 ! 18 ! 19 ! 20 ! 21 in several re- 
cent publications. It has been showniS that T c has a 
non-monotonic dependence on the thickness df of the 
ferromagnetic layers that provide information about the 
strength of the ferromagnetic exchange field and about 
the transparencies of the S-F interfaces. Approximate 
analytic formulas for T c have been derived for several 
limiting cases^ e.g. for thin or thick film thicknesses 
or for low or high interface resistances. Recently, 
Fominov et al. developed a numerical method to 
compute T c for diffusive S-F bilayers^ and symmetric 
F-S-F trilayers 1 ^ for arbitrary model parameters such 
as layer thicknesses and interface resistances. Such an 
approach is valuable when theory and experiments arc 
compared in detail with the aim to extract parameters 
as e.g. the ferromagnetic exchange field or the boundary 
transparencies. 

The possibilitjiSiiSii of a 7r-state (characterized by a 
stable phase difference of 7r between the superconducting 
order parameters) is now well established experimentally 
in S-F hybrid structures involving several superconduct- 



ing layers. Transitions between the and n states have 
been revealed in S-F-S junctions by the oscillations of the 
critical current when varying the temperature2£i2L2& or 
the ferromagnetic thicknessi 2 ^ ^ The transitions from 
the to 7r state may also be revealed^ by the presence of 
cusps in the dependence of T c on df. Because the cusps 
may be confounded with the oscillations of T c (df) them- 
selves, such a feature in the dependence of T c has been 
identified experimentally only recently^LS 

The presence of several ferromagnetic layers introduces 
a new degree of freedom, the relative orientation angle, 
9, between the magnetizations. The influence of the ori- 
entation on T c has been first studied theoretically in F- 
S-F trilayers in the Refs. and (these authors only 
considered parallel or antiparallcl orientations). The cal- 
culations for an arbitrary orientation were performed in 
Ref. 0. A dependence of the critical current oscilla- 
tions on the magnetization orientation has been also es- 
tablished theoretically in S-F-F'-S junctiona 33 i 34 i 35 i 36 and 
multilayered S-F junctions^ In Ref. |3 and|37] a switch 
between the and 7r states has been found from cal- 
culations of the Josephson critical current by changing 
the mutual orientation between the moments. The de- 
pendence of T c on the moment orientation (parallel or 
antiparallel) of t rilay ers has been studied experimentally 
in Refs. Ilall9ll20ll2ll A dependence on the domain state 
of the ferromagnet in a S-F bilayer was found in Ref. I22I 

Motivated by the recent experimental studies, we have 
developed a fast and effective method that is particularly 
suited for the numerical calculation of T c in diffusive hy- 
brid structures. An important part of this paper is to 
present details of this method and discuss the calcula- 
tions leading to some of the results presented in Ref. I32T 
Our method can be considered as a development of the 
method of Fominov et al., who in RefslSulCl have pre- 
sented calculations of T c of S-F bilayers and symmetric 
F-S-F trilayers with non-collinear magnetizations. We 
extend the calculations to the more general case of asym- 
metric trilayers in connection with the geometry consid- 
ered typically in experiments. Within our model, we also 
treat symmetric pentalayers, including the possibility of 
a phase difference 7r between the two superconductors. 
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This structure was recently studied experimentally in 
Ref. l32l From our T c calculations, we predict a switch- 
ing between 0- and 7r-junctions by the orientation of the 
ferromagnetic exchange fields in pentalayers consisting of 
a central Josephson junction, two superconductors sep- 
arated by a ferromagnet, sandwiched between two outer 
ferromagnets with exchange fields rotated relative to the 
central ferromagnetic layer. This kind of structure could 
be realized e.g. by fixing the moments of the outer layers, 
while rotating the moment of the central layer with an 
external magnetic field. 

The outline of the paper is as follows. In Section ITU we 
present the model of the F-S-F trilayer and the F-S-F-S- 
F pentalayer structures and outline the method that we 
use to compute the order parameter profile and T c of the 
structures. In Sections IIIII and llVl we present the results 
for the trilayer and pentalayer respectively. In Section 
[V]wc discuss some details of our numerical method. We 
summarize our work in Sect ion IVT1 Some of the technical 
details have been collected in the appendices. 



II. MODEL AND METHOD 

We shall restrict our considerations to diffusive hybrid 
structures and to temperatures T near the critical tem- 
perature T c . We employ a Green's function method in 
the quasiclassical approximation. The central quantity 
in this framework is the 2x2 spin-matrix anomalous 
Green's function /, describing superconducting correla- 
tions in the structure. The spin degree of freedom has 
to be kept due to the fact that the ferromagnets in prox- 
imity with the superconductors break spin rotational in- 
variancc. Thus, both spin singlet and spin triplet prox- 
imity pair amplitudes are present in the ferromagnet. We 
use a notation where the spin structure is described as 
/ = (f s + cr ■ f t )i<j y , where <r = (a x , a y , <t z ) are the three 
Pauli matrices. The pair amplitudes are the spin singlet 
component f s and the three spin triplet components de- 
scribed by the vector f t . The ferromagnetic regions are 
characterized by an exchange field with a fixed direction. 
In the case of rapid changes on the scale of the coherence 
length of the direction of the exchange field j^^MS or 
spin-active interface scattering^ long-range equal-spin 
triplet correlations are also induced. We refer to our re- 
cent paper92L224iMi and a recent review^ and references 
therein for a deeper discussion of the origin of these cor- 
relations. 

For diffusive structures the Green's function is 
isotropic to lowest order in l/pj£, where pf is the Fermi 
momentum and £ is the mean free path. Furthermore, for 
temperatures near the critical temperature the supercon- 
ducting gap is small A <C T c and the usual Green's func- 
tion is approximately equal to the normal state Green's 
function g ss — i7rsgn(e„), while the anomalous Green's 
function / is small, of the order of A. The relevant start- 
ing point in this case is Usadel's diffusion equation^ lin- 
earized for small A. We assume for simplicity that the 



spatial dependence in the structure is only along the in- 
terface normal, taken to be along the x-axis, see Figs. Q] 
and|SJ Then, the linearized Usadel equations have the 
formic 

(£>C - 2|£„|) fs = -2 7 rA + 2 l sgn( £n )J-f t (1) 
{Ddl x -2\e n \)i t = 2 lS gn(e„)J/ s , (2) 

where sgn(e„) is the sign of the Matsubara frequency 
e n = TtT(2n + 1), and we have used the short hand no- 
tation / = f(s n ,x). We assume that the exchange field 
J = 3(x) is non-zero in the ferromagnetic regions, while 
A = A(x) is non-zero in the superconducting regions. 
Each layer in the structure can have a different diffu- 
sion constant D. Note that we assume that the exchange 
field is small compared to the Fermi energy, J <C e/, in 
which case the quasiclassical theory can be straightfor- 
wardly applied. For strong exchange fields, a separate 
calculation has to be made4i 

The diffusion equation is supplemented with boundary 
conditions at each interface and at the outer surfaces of 
the structure. The boundary condition connecting the 
Green's function at xs on the superconductor side of the 
interface with the Green's function at xp on the ferro- 
magnet side of the interface is of the form first derived 
by Kupriyanov and Lukiche"\^i 

lipf'M = dsf(x s ), (3) 
Tbfrf'M = ±[f(x s )-f(x F )}, (4) 

where £ = yf D /2ttT c q is the coherence length and the 
parameters 7 and 7^ characterize the conductivity mis- 
match between the two sides and the boundary resis- 
tance, respectively. The sign in Eq. Q is positive (neg- 
ative) for a F/S (S/F) interface [for which the super- 
conductor occupies the space to the right (left) of the 
barrier] . Note that we use the prime as a short-hand no- 
tation for spatial derivatives at a certain point in space, 
e.g. f'(xs) = d x f{x)\ x=xs . At the outer surfaces of the 
structure, we require that the current trough the bound- 
ary must vanish, i.e. d x f = 0. 

Since the exchange field and the superconducting or- 
der parameter are spatially separated, we see that in the 
superconducting region Eqs. are decoupled. The 

triplet part J2J can be solved analytically, while the sin- 
glet part has a source term containing the order pa- 
rameter that satisfies the self-consistency equation 

A(x) ln-L = T £ (/.(*»,*) - ■ (5) 

T c0 ^ \ \e n \ ) 

In the ferromagnetic regions Eqs. Q-© arc coupled 
but the superconducting order parameter is absent and 
both equations can be solved analytically, which is de- 
scribed in detail in Appendix and [B] for the trilayer 
and pentalayer cases. The presence of the ferromagnetic 
regions arc in the process reduced to an effective bound- 
ary condition for the calculation of the singlet compo- 
nent in the superconducting region, which we confine to 
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< x < d s in the present discussion, as in Fig. ^ The 
boundary condition can in the general case be written in 
the form 



f'.(0) \_ kW ( /a(0) 



fs(d s ) 



(6) 



where k s = y / 2e n /D s and W is a 2 x 2 matrix. The 
non-locality of the boundary condition ©, i.e. the cou- 
pling of the two interfaces at and at d s , is a result of 
the coupling of the singlet and triplet anomalous Green's 
functions f s and ft in the original boundary conditions 
Eqs. ©-0 and the coupling of the diffusion equations 
for the singlet and triplet components in the ferromag- 
net by the exchange field. The matrix W depends on 
the Matsubara frequency, the parameters of the adjacent 
layers (thicknesses, exchange fields, diffusion constants), 
and the interface parameters 7 and 7^. The expressions 
for the components of W are derived in Appendix 1X1 and 
El for the trilaycr and pcntalaycr structures. The follow- 
ing method for calculating T c is however applicable for 
any matrix W, as long as the boundary condition for the 
singlet Green's function f s is of the form JJjJ • 

Consider Eq. JQ| in the superconducting region, i.e. 
for < x < d s where J = 0. By linear superposition we 
have& 

f s (e n ,x)='K G(e n ,x,y)A(y)dy, (7) 
Jo 

where the function G(e n ,x,y) is the solution of the dif- 
ferential equation 



\e n \) G(e n ,x,y) = -S(x -y), (8) 



2 xx 



subject to the boundary conditions © with f s (e n ,x) re- 
placed by G(e n , x, y). The solution of Eq. © is presented 
in Appendix IU1 With the help of the function G, the gap 
equation can be written as 



27rTE e „>o/o G(e n ,x,y)A(y)dy 



27TTV >o 



A(x), (9) 



where we used that the singlet Green's function f s {sn, x) 
[and therefore also G(e n ,x,y)] is an even function of e n . 
We see that one waj*& of solving the problem at hand is 
to discretize the spatial coordinate (x — ► Xk, k = 1...N) 
and find the critical temperature T c as the highest tem- 
perature for which the eigenvalue of the N x N matrix on 
the left hand side of Eq. equals one. The correspond- 
ing eigenvector gives the profile of the order parameter, 
A(x k ). 

There are several disadvantages of the method de- 
scribed above, all connected with the discretization of 
the spatial coordinate axis. In particular, it is cumber- 
some to reach acceptable numerical accuracy when T c is 
computed. We shall discuss these problems in detail in 
Section |V| 



Because of these draw backs, we develop a Fourier se- 
ries method that avoids the discretization of the spatial 
coordinate. The superconducting order parameter A(x) 
exists in the range < x < d s . We extend its domain of 
definition to the full real axis by adding an even-parity 
property and 2d s periodicity. Then, A(x) can be ex- 
panded in a Fourier series 



OO / s 



p=0 



(10) 



where the coefficients A p are defined as 



A„ 



2 -6. 



pO 



/ pTTx\ 



A(a;)cos dx. (11) 

V d s ) 

We show in Appendix [D] how to obtain an analytic ex- 
pression for the singlet amplitude f s in terms of the 
Fourier coefficients A p . Consequently, the gap equation 
can be written in the space of Fourier coefficients as 
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^rri/pAp = 0, 

p=0 

for integer I > 0, and where m; p are given in Eqs. i|D3|) - 
(|D4() . We solve the problem at hand by introducing a cut- 
off p c for the number of harmonics and find the critical 
temperature T c as the highest temperature for which the 
eigenvalue of the p c x p c matrix on the left hand side 
of Eq. l|12f) equals zero. The corresponding eigenvector 
gives the profile of the order parameter A(x) through the 
sum in Eq. i|10f> . 

In the following two section we use this method to 
compute T c for the S-F trilayer and pentalayer structures. 
In Section [V] we discuss the advantages of our method 
[Eq. 1(12)1] and compare with the other method [Eq. ©]. 



III. TRILAYER 

Consider the trilayer structure shown in Fig. ^ 
We study in this section the superconducting transi- 
tion temperature of such a trilayer. Our studies are 
motivated by the recent experiments on S-F layered 
struct ures) 12 i 13 i 14 i 15 i 16 i 17 i 18 i 19 i 20 i 21 including in particular 
the experiments in Ref . |3^ on the critical temperature of 
asymmetric F1-S-F2 trilayers. The theory fits of T c of 
the trilayers in Ref. |13 were obtained with the theory 
presented in the present paper. 

In the left ferromagnetic layer (Fi), the exchange field 
is aligned with the z-axis, while in F2 it lies in the yz 
plane and forms an angle 9 with respect to the z-axis. 
The origin of the coordinate system is taken at the Fi/S 
interface. The two layers Fi and F2 are characterized 
by their thicknesses (d/i, d/2), exchange fields (Ji, J2) 
and diffusion constants (Dfi, -D/2), while the supercon- 
ducting layer is characterized by its thickness (d s ), pair- 
ing interaction strength (i.e. the bulk material super- 
conducting critical temperature T c o), and diffusion con- 
stant (D s ). The diffusion constants are converted into 



FIG. 1: Geometry of the asymmetric F1-S-F2 structure. The 
moments Ji (in Fi) and J2 (in F2) may have different ampli- 
tudes and point in different directions (the relative orientation 
angle is denoted 9). 

coherence lengths £ = \J D /2-kT c q and we shall use the 
coherence length in the superconductor £ s as length scale 
in the problem. The Fi/S and S/F2 interfaces are char- 
acterized by the conductivity mismatches (71, 72) and 
interface resistances (761, 762)- 

The Usadel equations CO - © ar e solved as described 
in Appendix ^ to give the effective boundary condition 
matrix W for the trilayer. The matrix mi p of Eq. I|12(l 
is then given in terms of the elements of W as shown in 
Appendix [Dj 



A. Results 

In Fig. 12131 we present the influence of an exchange 
field on T c for an asymmetric F1-S-F2 trilayer (with 
dfi ^ df2). In the normal metal case (obtained by set- 
ting J = 0), the critical temperature is monotonically 
suppressed as the layer thickness dfi is increased, see 
Fig. [3 In the case of a ferromagnet, the exchange field 
induces an additional oscillatory behavior, closely con- 
nected to the spin mixing between up and down spins. 
As a result, T c is suppressed in a non-monotonic way. 
For a strong enough exchange field, the oscillation is so 
strong that superconductivity is suppressed at a critical 
thickness but can reappear at a larger thickness. This 
kind of non-monotonic dependence of T c was thoroughly 
studied 4 ! 5 ! 6 ! 8 ; 10 ! 11 ! 23 for F-S bilayers and symmetric F-S- 
F trilaycrs. The crossing of the two T c curves in Fig. |3| 
is due to the non-monotonic dfi dependence shown in 
Fig. 13 

The exact point where superconductivity disappears 
depends on many other parameters in addition to the 
strength of the exchange field. For example, the influence 
of the second ferromagnet 's thickness df2 can be under- 
stood in terms of the initial T c suppression at dfi = 0, 
which corresponds to the bilayer case. For sufficiently 



FIG. 2: Critical temperature T c of a trilayer versus the thick- 
ness of the left ferromagnet for several strengths of the ex- 
change field ranging from Ji = J2 = J = (the normal 
metal case) to strong J = 20T c o. The other layer thicknesses 
are d 3 = 2£s and d/2 = 0.5£s, the interface parameters are 
71 = 72 = 0.3 and 751 = 7^2 = 0.7, the exchange fields 
are parallel (9 = 0), and the diffusion constants are equal, 
D f i = D f2 = D s . 




J /T 

FIG. 3: Critical temperature T c of a trilayer versus the 
strength of the exchange fields Ji = J2 = J for a few layer 
thicknesses dfi of one of the ferromagnetic layers. The other 
parameters are the same as in Fig. 



large dfi the initial suppression is large enough that the 
subsequent oscillations with increasing dfi lead to dis- 
appearing and reappearing superconductivity, see Fig.0] 
Naturally, the initial T c suppression at dfi = is a non- 
monotonic function of df2, in analogy to the T c ((in- 
dependence. 

Another important parameter for the size of the T c 
variations, is the interface resistance. As seen in Fig. [5] 
superconductivity is suppressed in trilayers with good 
contacts (small 7&) for the model parameters chosen here. 
For example it is not possible to consider d s <C £, s , for 
which simplified calculations with a constant A can be 
made, and simultaneously consider good contacts 7^ — > 
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FIG. 4: Critical temperature T c of a trilayer versus the thick- 
ness dfi of the left ferromagnet for several thicknesses d/2 of 
the right ferromagnet. The exchange field is Ji = J 2 — 20T c n 
and the other parameters are the same as in Fig. [5] 
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FIG. 6: Critical temperature T c of a trilayer versus the misori- 
entation angle 9 between the exchange fields in ferromagnets 
Fi and F2. The exchange field is J 1 = J2 = 20T c n and the 
other parameters are the same as in Fig. [5] 
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FIG. 5: Critical temperature T c of a trilayer versus the inter- 
face resistance parameter 75 = 751 = 7^2 at a few thicknesses 
dfj corresponding to points on the J = 20T c o curve in Fig. |2] 
to the left, inside, and to the right of the T c — region. The 
other parameters are the same as in Fig. [5] 



FIG. 7: (a) The spatial dependence of the order parameter 
for several interface resistance parameters 7(, = 7(,i = 7^2 
on the solid curve (d/i = 0.1£s) in Fig. [S] (b) The Fourier 
components in Eq. I1UII . Note that A is normalized to the 
first component An, which remains unknown in a linearized 
theory. 



for reasonable conductivity mismatches (here 7 = 0.3). 
Certainly, for larger conductivity mismatches (smaller 7), 
T c is not suppressed as much and a smaller 7^ can be 
used. However, it is always important to keep in mind 
that T c is suppressed to zero in quite a large parameter 
space, including small d s and small 7b. 

In Fig. HJ1 we show the influence of the relative direc- 
tion of the exchange fields in the two ferromagnetic lay- 
ers. The dependence is monotonic, with the parallel ori- 
entation being the most destructive. We note (see also 
Ref. Hot) that for parallel or antiparallel exchange field 
orientations triplet correlations with zero spin projection 
on the local exchange field are present in the structure, 
while for intermediate orientations triplet correlations 
with non-zero projection are also induced. In order to 
describe the ^-dependence correctly, it is therefore im- 



portant to include ft, see Appendix [XI 

In Fig. |7{ a ) wc present order parameter profiles for 
four different values of 7^ on the solid line (dfi = 0.1£s) 
in Fig. El For a good contact (small resistance 7b) , the 
pair breaking becomes quite severe. The suppression is 
reflected as a growth of the Fourier components p > 1, 
see Fig.EJb). 



IV. PENTALAYER 

Consider the pentalayer shown in Fig. EI Experimental 
results for the critical temperature, including signatures 
of a transition from a zero- to a 7r-junction as function 
of the thickness of the central F-layer, were recently pre- 
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FIG. 8: The F2/S/F1/S/F2 pentalayer structure. We con- 
sider two types of misalignment of the outer exchange fields 
relative to the exchange field in the center layer. First, as 
shown here, J 2 is rotated by the same angle 6. The second 
possibility is when J2 is rotated in opposite directions, — 6 in 
the left F 2 and +6 in the right F 2 . 



sented for this structure in Ref. y2[ The theory fits of 
T c of the pentalayers in Ref. 112 were obtained with the 
theory presented in the present paper. 

The superconducting layers are considered geometri- 
cally identical with identical bulk material critical tem- 
peratures T c0 . In the central ferromagnetic layer (Fi), 
the exchange field is aligned with the z-axis, while in the 
right and left layers (F2) it forms an angle 9 with re- 
spect to the z-axis. We characterize the different layers 
by their thicknesses, exchange fields, and diffusion con- 
stants, with the constraint that the pentalayer should 
have certain symmetries with respect to the midpoint, 
see below. The present pentalayer problem can then be 
reduced to a trilayer problem with a new effective bound- 
ary condition at a fictitious outer surface at the center 
(x = 0). The two superconducting order parameters in 
the left and right S layers may differ in phase, which is 
reflected in the effective boundary condition. 

We shall consider two types of misorientation of the 
exchange fields in the outer layers relative to the center 
layer: the exchange fields J2 are rotated by +9 as in Fig. [3] 
(rotation type 1, +9/ + 9), or rotations by —9 and +9 in 
the left and right outer layers respectively (rotation type 
2, -9/ + 9). 

For rotation type 1 (+9/ + 9), when the phase differ- 
ence vanishes (0-junction case), the singlet component f s 
is an even function of x. Considering the parity of the 
exchange field J (J z — > J z and J y — ► J y ) and the Eqs. 
©-(0), we deduce that the ft z and f ty components have 
the same even parity. Thus, we impose the conditions 

(+9, +9), - jet : f' s (Q) = fl z (0) = 4(0) = 0. (13) 

On the other hand, when the phase difference is tt (tt- 
junction case), f s , f tz and f ty are odd functions of x and 
we impose the conditions 

(+9, +9), 7r-jct : f s (0) = f tz (0) = 4(0) = 0. (14) 

For rotation type 2 {—9/ + 9), the exchange field com- 
ponent J y is instead odd under x — > —x. For the 0- 




FIG. 9: Critical temperature T c of a pentalayer versus the 
center ferromagnet layer thickness 2d/i for several barrier 
transparencies. The curves come in pairs, solid line for the 
0-junction and dashed line for the 7r-junction, from top to 
bottom for 71,1 = ~/b2 = 7i> = {2, 1,0.8,0.7} respectively. The 
other parameters are ds = 2£s, d/2 = 0.5£s, Ji = J2 = 20T c o, 
6 — 0, 71 = 72 = 0.3, and D fl = D f2 = In- 



junction case, it implies that f ty is odd, while the other 
components are even just as above. For the 7r-junction 
case the parities are interchanged. The effective bound- 
ary conditions are 



{-9/ + 9), 0-jct 
{-9j + 9), 7r-jct 



f s (o) = &(o) = 0, 

fty(0) = 0, 



(15) 



f s (0) = f u (0) = 0, rifi , 

/ t ;(o) = 0. (lb > 



As shown in Appendix El the different boundary con- 
ditions yield different matrices W for the effective bound- 
ary condition (jSJl. 



A. Results 

The dependence of T c on the various parameters in 
the model is similar for the trilayer with left ferromag- 
netic layer thickness dfi and for the pentalayer with a 
phase difference between the two superconductors and 
a central ferromagnet layer thickness 2c?/ 1. The critical 
temperature is in fact equal for rotation type I, since 
the boundary condition at the center of the pentalayer, 
Eq. JT2J, is the same as for the outer surface of the tri- 
layer. Note, however, that the boundary condition for 
one of the triplets is different for rotation type II, sec 
Eq. (|15fl . The new ingredient in the pentalayer case is 
the possibility of a phase difference tt between the su- 
perconductors. The 7r-state can in simplified terms be 
understood as being due to the oscillatory behavior of 
the Green's function f s (s ni x) inside the central ferro- 
magnetic layer Fi. In an experiment, T c is given by the 
largest T c for each thickness and there will be a sudden 
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FIG. 10: (a) The order parameter profile for the parameters 
in Fig. El for j b = 0.8 at 2d/i = 0.4£s, for which the crit- 
ical temperature for the 0- and 7r-junctions are respectively 
T c « 0.3T c0 and T c = 0.16T c0 . (b) The corresponding Fourier 
components. 




FIG. 11: (a) The same curves as in Fig. [5] for 7 6 = 0.8. For 
thicknesses d/i between the two vertical lines the — > tt tran- 
sition can be tuned by the relative orientation of the exchange 
fields in Fi and F2. (b) The switch — » n appears at a critical 
angle 6 C ~ 84° for the thickness 2d/i = 0.65£s. For the larger 
thickness 2d/i = £s, outside the window indicated in (a), the 
largest T c is obtained for the 7r-junction. 



almost kink-like change in T c at the — > tt transition. For 
large oscillations, the transition becomes sharper. This 
is illustrated by changing the interface resistance % in 
Fig. [HI For good contacts and strong exchange fields, su- 
perconductivity can be destroyed at some critical thick- 
ness and then reappear at a larger thickness, just as in the 
trilayer case. For the pentalaycr, however, the 7r-phasc 
can pre-empt the 0-phase and superconductivity appears 
earlier compared to the trilayer as df\ is increased, see 
the curves for jb = 0.7 in Fig. 

An example of the order parameter suppression is 
shown in Fig. HUf ah corresponding to 7;, = 0.8 and 
2dfi = 0.4£s in Fig. [!|] The suppression of A at the 




(b) 

"0"-phase 



J=10T, 



"7i"-phase 



FIG. 12: (a) Phase diagram of the — > tt transition in the 
window indicated in Fig. lUf a). In (b) we show the phase 
diagram for a smaller exchange held, J = 10T c o. The range 
of thicknesses d/i for which there is a switching by changing 
is larger in this case. 



interfaces is more severe for phase difference tt and the 
0-junction is stabilized, i.e. has the largest T c as seen in 
Fig.H 

In the region close to the — > tt transition, it is possible 
to switch between the 0- and 7r-phases by changing the 
relative orientation of the exchange fields. We note that 
this possibility was already deduced from calculations of 
the Josephson critical current in the papers an d 
considering different geometries. We illustrate this effect 
in Fig. 1111 switching is possible in between the vertical 
lines in Fig. lUf a). Since experimentally, T c is given by 
the largest T c for each 9, the — > tt switch would show 
up as a sudden almost kink-like change in T c with the 
variation of 9, as shown in Fig. Illf bh We present in 
Fig. El the phase diagram of the junction in the region 
around the window indicated in Fig. Illf a). The window 
inside which a — > tt phase change can be induced by 
the orientation angle 9 is larger for a smaller exchange 
field since the T c -oscillation period is longer in this case. 
We see this effect by comparing the J = 20T c o case in 
Fig. lUf a) to the J = 10T c o case shown in (b). 

It has been found&2iS for the bilayer and trilayer cases 
that T c can become a multiple valued function of e.g. the 
thickness of the ferromagnct. We show this type of be- 
havior for the pentalayer case in Fig. 1131 (upper panel). 
The non-monotonic dependence of T c is similar to the 
case of a clean thin film in an in-plane magnetic field^4S 
and to thin films of superfluid 3 He |50ll5lT |. For these 
clean systems it has been proposed that an inhomogc- 
neous superconducting state can be formed. In a dirty 
system, such inhomogencity seems very unlikely and it 
has instead been proposed that the back-bend signals 
the possibility of a first-order transition in the system^ 
First-order transitions are however beyond the scope of 
the present paper. Instead we point out that for the pen- 
talayer case, the 7r-phase becomes favorable in the same 
region of thicknesses as where there is a back-bend for the 
0-junction. The back-bend behavior for the 0-junction, 
and the interfering 7r-phase, occurs also as function of the 
exchange field misoricntation angle 9, sec the lower panel 
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FIG. 13: Upper panel: In the region where T c is suppressed 
to zero, the curve T c (d/i) can contain a back-bend. This 
latter could signal the occurrence of a first-order transition, 
which is however beyond the scope of the present theory. For 
the pentalayer, the 7r-phase can interfere and the first order 
transition might be avoided. The parameters are d a = 2£ s , 
d f2 = Q.2i s , 71 = 72 = 0.35, 7w = j b2 = 0.4, Ji = J 2 = 10T c0 
and 9 = 0. In the lower panel we study the dependence on 
the exchange field for two particular thicknesses dfi in the 
upper panel. Clearly, the back-bend behavior can occur also 
as function of the exchange field misorientation angle 9. 

in Fig. 1131 Interestingly, there is a discontinuous drop in 
T c at the — ► 7r transition when 9 is tuned from around 
20° down to 10°, see the solid and dashed lines in Fig. 1131 
For very large thicknesses dt\ the predominant su- 
perconducting correlations that penetrate Fi and con- 
nect the two superconductors are the long-range non- 
oscillatory triplet components of ft. As a consequence, at 
large dfi, T c becomes a monotonic function of dfi. The 
difference in T c between the 0- and 7r-phases is, however, 
quite small, see Fig. ^] The junction is stabilized at 
large dfi either at or at ir phase difference, depending 
on the way the exchange fields of the two outer ferro- 
magnetic layers are rotated relative to the center layer (a 
similar effect associated with the chirality of the rotation 




FIG. 14: For thick center films (large dfi) the communication 
between the two superconductors is taken over by long-range 
non-oscillatory equal spin triplet correlations and the junc- 
tion is stabilized at phase difference or tv depending on the 
exchange field orientation: the upper panel (7r-junction for 
9 = 90°) is obtained for exchange field rotation type 1 (+9 in 
both the left and right outer ferromagnets F2, as illustrated 
in Fig. IS) , while the lower panel (0-junction for 9 = 90°) is 
obtained for rotation type 2 (—9 in the left F2 and +9 in the 
right F2). The difference in T c between the and n cases 
is quite small for large dfi. Upper inset: the differences in 
T c for various exchange field orientations (solid lines 9 = 0, 
dashed lines 9 = 90°) are due to the interaction between the 
ferromagnetic layers Fi and F 2 . Lower inset: spatial depen- 
dence inside the central 6£ s thick Fi layer of the long-range 
triplet ft y induced for 9 = 90°. The parameters are d„ = 2£ s , 
d f2 = 0.5£ 3 , dfi = 3£ s , 71 = 72 = 0.3, 7bi = 7 62 = 0.8 and 
Ji = Ji = 10T c0 . 



has been found in Ref . |3]] from calculations of the critical 
current in S-F multilayered junctions). Wc consider two 
types of rotation: the exchange fields in the outer ferro- 
magnets are rotated by +9 as in Fig. |S] (rotation type 1, 
+9/ + 9), or rotations by —9 and +9 in the left and right 
outer layers respectively (rotation type 2, —9/ + 9). The 
only difference between the two phase differences and tt 
is the parity property of one of the triplets: f ty is an even 
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(odd) function of x for the 0-junction (7r-junction) for ro- 
tation type 1 (+0/ + 0). The parity properties of f ty for 
and 7r phase differences arc reversed for rotation type 
2 (—0/ + 6). When f ty has odd parity it is smaller com- 
pared to the even parity case, which leads to a smaller 
suppression of the singlet f s , i.e. less pair breaking, and a 
higher T c . We therefore have a 7r-j unction at large dfi for 
rotation type 1 (upper panel in Fig. ll4f) . and a 0-junction 
for rotation type 2 (lower panel in Fig. ITU). We show the 
spatial dependence of the long range (in the central F± 
layer) triplet Green's function 
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in the lower inset of Fig. PHI 

Experimentally, the transition from — *• 7r was stud- 
ied until now by varying the thicknesaii22i2fli2ii22 of 
the ferromagnet in S-F-S junctions, or by varying the 
temperature^SSiSLS which is more practical since the 
transition is seen in the same device. Here we have 
studied another possibility to switch from the to the 
7r state within the same device, namely by continuously 
changing the relative orientation of the ferromagnetic 
moments. Our results are qualitatively consistent with 
the results obtained within Josephson critical current 
calculations The feasibility of controlling the orien- 
tation of the moments has been proven experimentally 
through the investigation of F-S-F trilayers for different 
moment oricntationsii£ii2i22i2i 
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V. DISCUSSION OF THE NUMERICS 

Let us discuss some delicate problems that need to be 
addressed when T c is computed in inhomogeneous struc- 
tures. In particular, we will compare the two methods of 
computing T c : Eq. (|12|l which we call the Fourier method 
and Eq. © which we call the grid method. 

The most important problem to address in any cal- 
culation using Usadel's approximation is the fact that 
the Matsubara sum in Eq. (jSJ) is intrinsically slowly 
convergent, as compared to calculations done with the 
more general Eilcnberger approach. As we show in Ap- 
pendix^ the difference / s (e n )— 7rA/|£r ra | appearing in the 
gap equation (JSJ is at high-energies proportional to 1/ e\ 
for inhomogeneous systems. This can be contrasted with 
an Eilcnberger approach, where the high-energy asymp- 
totic is 1 / l^n | 3 - It is therefore always necessary to ex- 
tend the Matsubara sum to high energies when the Us- 
adel approximation is employed, see the dashed line in 
Fig. El In the example we need a technical cut-off of 
order 1000T c o to compute T c with an accuracy of 1%. 
However, since the high-energy form of / s (£n) is known 
(see Appendix [EJ it is in principle possible to circumvent 
the problem by treating the high-energy tail separately 
and sum the Matsubara sum to infinity. We have done 
that within the Fourier series approach, see Appendix IP! 



FIG. 15: Upper panel: critical temperature versus the techni- 
cal cut-off e c . To achieve good accuracy for the critical tem- 
perature we need e c of order 1000T c o (dashed line). When 
the high-energy tail is summed to infinity, as described for 
the Fourier method in Eq. ijFlOfl . the convergence is more ac- 
ceptable (solid line). Lower panel: the eigenvalue of the gap 
equation 11211 versus temperature for two different thicknesses. 
The zero-crossing determines T c . When T c is suppressed, A(T) 
can become a flat function of T which makes it important to 
compute A with high accuracy to avoid numerical errors in T c . 
The parameters in the upper panel were chosen as in Fig. 1131 
0-junction, at 2d/i = 0.7 and 8 = 0. In the lower panel the 
two thicknesses are indicated in the legend. 



A more acceptable cut-off of order 100T C is then enough 
to achieve excellent accuracy, see the solid line in Fig. 1151 

There are several other factors that, together with the 
slow convergence of the Matsubara sum, conspire to make 
it non-trivial to achieve acceptable accuracy, especially 
when T c is small compared to T c q. The critical temper- 
ature is computed by finding the temperature for which 
the eigenvalue A of the gap equation is zero [Fourier 
method, Eq. l(T2)l] or one [grid method, Eq. ©]. The 
function A(T) can become a very flat function of T in the 
region where T c0 is small, see the lower panel of Fig. ITSl 
Any error made in the calculation of A can therefore be 
magnified to a larger error in T c and it becomes increas- 
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ingly critical to compute A with high accuracy as T c is 
suppressed. 

The above two technical problems arc particularly hard 
to circumvent within the grid method. First of all, the 
need to include high energies up to a technical cut-off 
e c imposes a condition on the grid spacing 5x. At high 
energies the function G(e n , x, y) is typically peaked in the 
region x ~ y 

G(e n » T c0 ,x,y) ~ Mf^ e -k.M\*-v\ (18) 

where k s (e n ) = \j2e n /D. It is therefore necessary to 
choose 

*£« I = ,/^° (19) 

to resolve this dependence. Since we need a cut-off 
around 1000T c o, because of the slow convergence within 
the Usadel approach, we need a grid spacing of order 
0.01£s or finer. The matrix in Eq. JjJJ must therefore typ- 
ically be of the order of a few hundred elements square, 
which severely slows down the numerics. 

One reason for the importance to resolve the peaked 
form of G(e n ,x,y) is due to the interchange in order of 
the Matsubara sum and the integration over y in Eq. 
We write Eq. © as 

[ * K(x,y)A(y)dy = A(x), (20) 
Jo 

and compute each element of the matrix K(x,y) by sum- 
ming over e„. The asymptotic form of the diagonal is 
however G(e n ,x,x) oc and the Matsubara sum is 

not convergent. This is in principle irrelevant for the cal- 
culation of T c because T c only depends on the eigenvalue 
of the matrix, which is a quantity given by the Matsub- 
ara sum integrated over y. Note that when Eq. l|18Jl is 
integrated over y, a factor l/k s appears in the primitive 
function of the exponential and the asymptotic form is 
1/|£„|, which is (by construction) cancelled by the sum 
over l/|evi| in the denominator of K(x,y), see Eq. 
Numerically, however, the integral over the discrctized 
coordinate y can only be computed with some accuracy 
given by the grid spacing Sx. The error made in comput- 
ing the integral is transferred into an error in the eigen- 
value A which, as described above, can result in an error 
in T c magnified by the flatness of the A(T)-dependence. 

To circumvent the problems described above, one must 
predict the high-energy tail to avoid cut-offs larger than 
~ 100r c o. Within the grid-method that means comput- 
ing the derivative of A.(x), i.e. to introduce an approxi- 
mate formula for the derivative on a discretized grid. But 
that also introduces numerical errors and the grid must 
still be dense, which means that the matrix K(x,y) re- 
mains large and the calculation with the grid method is 
always very slow and susceptible to numerical errors. 

All the problems related to the discretization of the 
spatial coordinate are avoided within the Fourier series 




FIG. 16: Critical temperature versus the number p c of in- 
cluded Fourier coefficients in Eq. ffllUB . The variation in abso- 
lute numbers are shown by circles (vertical scale to the left), 
while the variation in relation to the corresponding value for 
T c at a high cut-off p c — 100 is shown by squares (vertical scale 
to the right). The even-odd variation is due to the choice of 
parameters: the junction is almost symmetric and the even- 
number Fourier components corresponding to symmetric cos- 
functions contribute more to T c . The model parameters were 
chosen as in Fig. 1151 upper panel. 

approach, since G(e n , x, y) is analytically integrated over 
x and y in the course of the derivation of the matrix mi p 
in Eq. I|12(l. see Appendix iDl Moreover, the high-energy 
tail is easily predicted analytically, see Appendix IF1 It 
is typically sufficient to include only the 20 first Fourier 
components in the calculation of T c , see Fig. The 
matrix mi p is therefore small, the high-energy cut-off of 
the Matsubara sum can be chosen reasonably small, and 
very high accuracy is achieved while the speed of the 
calculation remains very high. 



VI. SUMMARY 

In conclusion, we have studied the change of the super- 
conducting critical temperature, T c , in asymmetric tri- 
layers Fi-S-F 2 and symmetric pentalayers F 2 -S-Fi-S-F 2 
with any relative orientation angle between the magneti- 
zations of Fi and F 2 . For both cases we have presented 
phase diagrams, showing T c as function of the misorien- 
tation angle, 9, and as a function of the ferromagnet layer 
thicknesses. We have investigated the interplay of long- 
range triplet components and Joscphson coupling in the 
pcntalayer geometry. We have demonstrated the possi- 
bility to switch between the and 7r states by controlling 
the relative orientation of the F moments in a pcnta- 
layer structure. This behavior may be appealing for the 
experimental study of the — > tt transition. We have 
presented details for a general method for the computa- 
tion of T c and the dependence of the order parameter 
on the spatial coordinates in diffusive hybrid structures. 
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With this technique, the accuracy as well as the speed 
of the numerics are immensely improved compared with 
previously used techniques. 
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APPENDIX A: DERIVATION OF W FOR THE 
TRILAYER 

In this Appendix we provide the details of the calcu- 
lations leading to the effective boundary condition (JSJ) 
obeyed by the singlet component in the superconduct- 
ing region in asymmetric F1-S-F2 trilayers with an ar- 
bitrary mutual orientation between the magnetizations 
in Fi and Fj. Except for this latter component, it is 



possible to derive analytically the spatial dependences of 
all the components of the anomalous Green's function / 
close to T c (next Section). In Section TA 21 we determine 
from the consideration of the boundary conditions @-(QJ 
the matrix W that enters the expression for the effective 
boundary condition (J5J. 



1. Spatial Dependences 

In the superconducting layer, the triplet vector f t obeys 
a homogeneous differential equation [Eq. J2J with J = 0] 
which is straightforwardly solved: 

ft = ccosh(/s s .T) + d sinh(fc s x) (Al) 

with c and d constants. 

For a fixed exchange field in each F layer, the sys- 
tem of coupled Eqs. Q-© can be easily solved in the 
ferromagnetic regions. After application of the bound- 
ary conditions at the outer surfaces, the solutions can be 
written in the form— 



f*) = a eCOsh[k el (x + dfi] \ J +a cosh[/c i(a; + rf/i)] ^ ? ^ (A2) 



for the Fi layer, and 



t s I = b e cosh \k e o(x — d s — dto)] ( , n <-^, ■ n -\ I + &o coshffcn2(^ — d 8 — dfz)] ( . n - [A3) 

i t J £ L 6 y 1 n y e(cos0z + sm#y) J ' 1 u v 1 n y smflz - cos6>y / ' 



k± q = J{2e n ±2iJ q )/D fq , (A4) 



e=± 

for the F2 layer. Here wc have defined 



k 0q = yJ2e n /D fq , (A5) 
with the index q = 1 or 2 referring to the Fi or F2 layer. 

2. Determination of W 

The constants aj and bj (j — ±,0), c and d are determined with the help of the boundary conditions ©-(01 
considered for the two S/F interfaces. Writing these conditions for the triplet components only, we have 

UKx Sq ) = l^fj'dxp,), (A6) 

fl{x S q) = fl(xFq) + VqlbqtfqflixFq), (A7) 

with I — t y ,t z . Note that 771 = +1 and 7/2 = — 1- Similarly, wc get for the singlet amplitude 



Zsf's(XSq) = lqtfqfs(XFq), 
fs(xSq) = fs(xFq) + Vq1bq€fqf' s {XFq)- 



(A8) 
(A9) 
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Here, xfi and x$i are the coordinates on the two sides of the Fi/S interface at x\ = 0, while xf2 and Xs2 refers to 
the S/F 2 interface at X2 = d s . From Eqs. (|A9() and (|A7() for the first interface (Fi/S), we obtain the system 

f s (xi) = ^a e A, (A10) 

e=± 

Cy = UqAq, (AH) 

c z = ^2sa s A e , (A12) 

e=± 

with the quantity 

Aj = cosh(kjidfi) + jbikji^fi sinh(kjidfi), (A13) 

where j = ±, 0. The matching of the different components with the conditions l|A9(l and l|A7(l yield at the second 
interface (S/F 2 ) 



f.(x 2 ) = 5> £ S £! (A14) 

e=± 

c y cosh(k s d s ) + d y smh(k s d s ) = eb £ B £ sin 9 — b a Bo cos9, (A15) 

c z cosh(k s d s ) + d z sinh(fc s d s ) = £& e £> £ cos 6* + & <8o sin0, (A16) 



with 



= cosh(/cj 2 d/2) + 7b2%2^/2 sinh(fej 2 d/2) (A17) 
defined in a similar way as the quantity Aj . Then, the boundary conditions l)A6jl yield the system 

d y = a C , (A18) 

d z = ^ea £ C £ , (A19) 

e=± 

for the Fi/S interface, and 

Cy sinh(fc s <i s ) + d y cosh(fc s d s ) = b^T) cos# — sb £ T> £ sin9, (A20) 

e=± 

c z sinh(fc s d s ) + d z cosh(fc s <i s ) = —b T> sin — sb £ V £ cos 0, (A21) 



for the S/F 2 interface, with 

= 7i fc ji£/i sinh(fc i i<i / i)/fc s ^, (A22) 
£>j = 72^2^/2 sinh(fc j2 d /2 )/fc s ^. (A23) 

The next step consists of eliminating the coefficients c y , c zi d yi d z , a , 60 from the former equations. We obtain the 
system 

£)ea e £ e = 5>& £ e £ , (A24) 

e e 



where 



£ £ = K (A £ -C £ )[cosh(k s d s ) -sinh(k s d s )} cos 9, (A26) 

T £ = K e (B - V ) sin 2 6> + K Q (B £ -V £ ) cos 2 0, (A27) 

£ e = A' £ (6o+^o)sin 2 + Xo(6 e -+^e)cos 2 0, (A28) 

H e = A' (A +C e ) [cosh(fc s 4) + sinh(fc s d s )] cos6>, (A29) 
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with 



Kj = [B-jCo + VjAo] cosh(k s d s ) + [BjAq + VjCo] sinh(Ms)- (A30) 
Compiling Eqs. IjAlOfl and (|A14|) with Eqs. (|A24|) - (|A25|) . we get the expressions for the amplitudes a e and b E 

(B+l-,- £ + BJI + ,- s ) f 8 (xi) + eA- e {T-Q+ - T+Q-) f s (x 2 ) 



with 



Finally, Eqs. 



be = 



J 

{A + l- e .- + A-I- E , + ) fs(x 2 ) + £g-e (£-H+ - E+H-) fsjxi) 

J 



T e E e i - Qe'He 



J = A+B + I-,-+A-B-l + , + + A + B-l + ,-+A-B+l-. 
yield the system 

Ufa) - -fc a &5> e X> c , 



(A31) 
(A32) 



(A33) 
(A34) 



(A35) 
(A36) 



which can be rewritten in the form 



f s {xi) \_ k (W n W 12 \f f s (xi) \ 
f'sfa) ) ~ s \ W 2 i W22 ) \ f s (x 2 ) J ' 



with 



W n 
W 22 
W 12 
W 21 



C+ {B+l- 



BJI + ,-)+C- (B+T-.+ + BJI +A 



J 

V+ (A+1-- + A-1-.+) + V- (A+1+.- + A-1+. 

J 



(^G+ - (A-C+ - A+C-) 

J 

(£-H+ - S+H-) (B-V+ - B+D-) 
J 

Using the expressions <|A26|l - <jA29|) . one can notice that in fact W\ 2 = —W 2 \ with 

2Kl cos 2 6 {B-V+ - B+V_) (A-C+ - A+C- 



W 12 = 



J 



(A37) 

(A38) 
(A39) 
(A40) 
(A41) 

(A42) 



For an asymmetric trilayer F1-S-F2, the diagonal coefficients Wn and W 22 of the matrix W differ in general. In the 
special case of a symmetric trilayer F1-S-F1, we have Wn = W 22 . 



APPENDIX B: DERIVATION OF W FOR THE 
PENTALAYER 



Due to the symmetry of the geometry, we need to de- 
termine the components of the anomalous Green function 
/ only in half of the pentalayer, e.g. in the domain x > 0. 
The problem is mapped back onto the asymmetrical Fi- 
S-F 2 trilayer problem previously considered in Appendix 
1X1 Because we have chosen a different origin for the sys- 
tem of coordinates, the Fi/S and S/F2 interfaces are now 



r 



located at the positions x\ — df \ and x 2 = d s + Due 
to the shift in coordinates, we have used the expressions 
(|A1|) in the S layer and (|A3|) in the F2 layer with x re- 
placed by x — dfn 

For the rotation type 1, the spatial dependences of the 
singlet and triplet components of / in the left Fi layer are 
in the 0-junction case the same as in Eq. (|A2|) after the 
shift of coordinate. In the 7r-junction case, the boundary 
conditions at the (fictitious) outer surface x = have 
changed, and the spatial dependences in Fi arc now given 
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by: 



where we introduced the notation 



e=± 



a £ sinh [k e \x] 



ez 



-a sinh [ko \x] 



The new boundary conditions (|13|1 and l|14(l at x = 
do not affect the definition of the former quantities Bj 
and Dj. On the other hand, changes occur in the def- 
inition of the quantities Aj and Cj (where j = e or 0). 
In the 0-junction case, the coefficients Aj and Cj remain 
unchanged, while in the 7r-junction case they are defined 



■4? 
c, 



sinh(fcjid/i) + jbikji^fi cosh(fejid/i), (Bl) 



7i%i£/i cosh(%id/i)/fc s ^ s 



(B2) 



For the rotation type 2, the components of / in Fi have 
a different spatial dependence as a result of the conditions 
(|15f) or (XBt - They are expressed as 



f * ^ = ^ cosh [ft e ix] ^ £ \ ^ +a sinh [fc ix] ^ ? 

in the 0-junction case, and 

■j? ^ = a e sinh [fc e ix] ^ g \ J +clq cosh [Aft 1&] ^ 



in the 7r-junction case. As for rotation type 1, changes 
occur in the definition of the quantities Aj and Cj (where 
j = e or 0) for the rotation type 2. In the 0-junction case, 
the coefficients A £ and C e have the same expression as in 
Appendix^ while Aq and Co are now given by 

A Q = sinh(fc i<i/i) +7bifcji0i cos h(fc id/i), (B3) 
Co = 7i fc oi£/i cosh(fcoid/i)/fc s £ s . (B4) 

In the 7r-junction case, the quantities Aq and Co are de- 
fined in the same way as in Appendix^] while A e and 
C e are written as 

A e = sinh(fc e id/i) + 7&ife e i£/i cosh(fe id/i), (B5) 
C £ = 7ifc E i^/icosh(A; E id/i)/fc s ^ s . (B6) 

Except for these modifications in the definition of the 
quantities A and C, the remaining calculations are ex- 
actly the same as in the asymmetric trilayer geometry 
and we can use the final expression derived in Appendix 
lAl for the matrix W in the symmetric pentalayer struc- 
ture. 



X\[x) = cosh(fc s x), 

X 2 (x) = sinh(fc s x), 

Yi (x) = cosh (k s [x - d s ] ) , 

Y 2 (x) = sinh (k s [x — d s ] ) . 



(C2) 
(C3) 
(C4) 
(C5) 



The coefficients L c , L s , R c and R s depend on the location 
y of the source term in Eq. (JSJ. The source is taken into 
account by the conditions 



G(x,y)\ x = y + = G(x,y)\ x=y - , 



(C6) 



and 



d x G{x,y)\ x=y+ - d x G(x,y)\ Q 



-k 2 Je n , (C7) 



where y + and y~ denote the limits y — » x from above and 
below, respectively. Eqs. (|C6|) - (IC7t give two relations 
between the coefficients in Eq. (IClf) . Two additional re- 
lations are provided by the boundary conditions at the 
edges of the superconductor, which read 



dxG(x,y)\ x=0 
dxG{x,y)\ x=ds 



KW 



G(0,y) 
G(d s ,y) 



(C8) 



These conditions are consistent with the boundary con- 
ditions © obeyed by the singlet amplitude f s (x). Com- 
piling Eqs. (|C1|) - (|C8|) . we obtain the coefficients 



L c (y) 
R c (y) 
L s (y) 
Rs(y) 

where 



k s 

£n£ 

ks 

£ n C 

k s 
e n C 

k s 
p r 

on**' 



[Yx(y) + W 22 Y 2 {y) - W 12 X 2 (y)} , 
\Xi{y) - W 21 Y 2 (y) + W n X 2 (y)] , 
WnY^y) + WiMy) + det{W)Y 2 (y) 
WMy) + W 22 X x {y) + det{W)X 2 (y) 



C = W 12 -W 21 + (W n -W 22 )cosh(k s d s ) 
+ 1 - det(W) smh(k s d s ). 

We note that the dependence on the Matsubara fre- 
quency e„ enters through k s and the four elements Wn, 
W 22 , W12 , and W 22 of the 2x2 matrix W in the boundary 
condition. 



APPENDIX C: DERIVATION OF G(e n ,x,y) 

In analogy with Ref . |H Eq. © is solved by making the 
following ansatz 

r(n , ,a _ / L c (y)Xi(x) +L s (y)X 2 (x), x<y, , , 
^ x ' y > - { R c (y)Y 1 (x) + R 8 (y)Y 2 (x), y < x, 



APPENDIX D: DERIVATION OF m ip IN EQ. JT3J) 

We insert the expansion (fTTTfl into Eq. J7J), use the ex- 
pression for G(e n ,x,y) derived in Appendix IU1 and per- 
form the integration over the spatial coordinate y. We 
obtain the singlet amplitude f s (s n ,x) in terms of the 
Fourier coefficients A B 
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f s (e n ,x) = — ^2 A P (3 P I C cos ( ^— j 



+ [VF21 + (-1) P W 22 ] Xi(a;) - [Wh + (-l)*Wi 2 ] Yx(x) 
+ det(W) l(~irx 2 (x) - Y 2 (x)] } , (Dl) 

where /3 p = 1/ [l + (7rp/fc s <i s ) 2 ] and the functions £, A^x), Yi(a;), and Y2(a0 were introduced in Appendix IU1 

We insert this expression in the gap equation @ and project in Fourier space, i.e. we multiply by cos(irlx / d s ) and 
integrate over x. As a result, we obtain a linear system for the Fourier components A p , with row / > given by 



+00 



2J m; P Ap = 0. 

p=0 

The off-diagonal elements (I 7^ p) have the form 

mi p = AttT ^2 — hp fliPp, 

e„>0 £n 



while the diagonal elements (Z = p) are given by 

mu = (1 + S l0 ) In — + 4ttT ^ — 



£„>0 



(D2) 



(D3) 



(D4) 



where 



[W11 - (-l) i+p ^ 22 + (-lyWla - (-1)^21] sinh(fc s 4) + det(W) {(-!)*■ + (-1)' - [1 + (-l) i+p ] cosh(M.)} 



1 — 



1M 



The relation W12 = —W21 between the off-diagonal ele- 
ments of W found in Appendix^implies that the matrix 
m is actually symmetric, i.e. mi p = m p i (see expressions 
(ID3(I and i|D5|0 . This property guarantees the existence 
of real solutions of the eigenproblcm Q12JI. 



where g is a 4 x 4 matrix in combined particle-hole 
and spin spaces, fj (j = 1, 2, 3) are the Pauli ma- 
trices in particle-hole space, and A is the gap function 
(A = (ia y )fiA if A is real). Eq. (|E 1 ft is supplemented 
with a normalization condition 



APPENDIX E: HIGH-ENERGY ASYMPTOTICS 

We present and compare the asymptotic high-energy 
behavior of the quasiclassical Green's function in the dif- 
fusive limit within the Usadel approximation to the more 
general case described by the Eilcnberger equation. Since 
the present discussion is independent of the presence or 
absence of a weak exchange field J< in the system, 
we leave it out. 



1. Diffusive Limit 

The Usadel equation^ for arbitrary temperatures (not 
necessarily close to T c as in the rest of the paper) is 



ie n T 3 - A, g 



+ -d x (gd x g) = 6, (El) 

7T 



-2 2i 

g = -it 1. 



(E2) 



Further details concerning the structure of the Green's 
function with the present notation can be found in Ref.l43l 
(see also Ref . l52|) . 

At high energies the order parameter and the deriva- 
tive term are small, 



A ~ T c0 <C e„, 
D/e ~ T c0 « e r , 

and we expand the Green's function 



9 



(0) 



,(D 



i(2) 



(E3) 
(E4) 



(E5) 



where the term g^ is of order (T c0 /e n ) k . To lowest order 



we have 



(0) 



-TT 2 !, 



(E6) 
(E7) 
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with the solution 



2. Arbitrary Mean Free Path 



g(°) = (-wr)sgn(e„)f 3 . 



In hrst order we obtain 



(i) 



0. 



(E8) 

(E9) 
(E10) 



Since g^ is proportional to f 3 , the second line can be 
used to move g^ 1 ' to one side of the commutator on the 
left-hand-side of the first line. Wc obtain 



(Ell) 



2^e„f3g (1) 

The solution is purely off-diagonal in particle-hole space 
»W = t^l (r 3 At 3 -&)= 'A. (E12) 



2i\e. , 

In second order we have 



7T 

2 



(E13) 
(E14) 



After a short calculation, similar to the calculation in 
first order, we obtain 



9 



(2) 



f 3 A 2 



(E15) 



Note, in particular, that there is an off-diagonal term 
proportional to \je 2 a for inhomogeneous systems. 

The off-diagonal part of the Green's function has ac- 
cording to the above the asymptotic form 



f{£n,x) 



irA(x) irDd 2 x A{x) 



-0 



\e n \ 



, (E16) 



which we now use to discuss the gap equation. The gap 
equation 



A(x) = XT 



E 

|£„|<W p 



/(e n , x), 



(E17) 



contains a log-divergency and it is necessary to introduce 
a cut-off uj p . But by the well-known procedure (see e.g. 
Ref. 153^ , the interaction strength A and the Matsubara 
sum cut-off tUp can both be eliminated by adding and sub- 
tracting the leading high-energy term in Eq. I|E16(I . The 
gap equation then has the form in Eq. ijjjj. The Matsub- 
ara sum converges, with a high-energy asymptotic tail 
oc \je 2 n according to Eq. (|E16|) . and can be extended to 
infinity. In practice, a technical cut-off e c is introduced 
that should, however, be high enough that the results of 
the calculation are cut-off independent. 



We compare the above results obtained within the Us- 
adel approximation with the corresponding high-energy 
behavior obtained within the Eilenberger approach. The 
Eilcnbcrgcr cquation^^^ reads 



is n h - A - a imp ,g 



iv F ■ Vg = 0, 



(E18) 



with impurity self energy <7i m p, and where vp is the Fermi 
velocity. The normalization condition g 2 = ~ir 2 l holds. 
We include non-magnetic impurity scattering within the 
self-consistent i-matrix approximation, for which the im- 
purity self energy is 



,(s) =ci(s,s), 



(E19) 



where c is the impurity concentration, and s is a param- 
eter that specifics the position of the momentum on the 
Fermi surface. The i-matrix is given as the solution of 
the equation 



t(s, s') = u(s, s') + 

'u(s,s")M F (s")g(s")i(s",s') 



(E20) 



where we have omitted for brevity all variables except 
the Fermi-momentum. Here, u(s,s') = w(s,s')l is the 
impurity scattering potential, and denotes a Fermi 

surface average over ,s' . 

We expand g as in Eq. I|E5(I . The zeroth order term for 
the Green function is given analogously to the discussion 
for the diffusive limit by 



*(0) _ 



= (-wr)sgn(e n )f 3 . 



(E21) 



For the higher orders we need to expand the impurity 
t-matrix in the parameter (T c o/s n ) 7 



i = m + + + 



(E22) 



and similarly for the impurity self energy. Introducing 
the operator 



D(s, s') = S(s - s')l - u(s, s')Af F {s')g 



(0) 



(E23) 



the t-matrix equation for the zeroth order term takes 
the form 



(E24) 



^( S ,s")i (( V',s')) n =u(s,s'). 
With the inverse operator defined by 

(i?- 1 ^, *")£(*".*')) h = s{s-s')i 

the formal solutions are given by 
P\a,a') = ^- 1 (s,s")u(s",s')) „ 



(E25) 



(E26) 



i (1 \s,s') = (^( S , S ")AA F ( S ")5 (1) (s")i (0) (s", S ')) s „ 
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From Eq. HE26f> wc obtain. 



Matsubara sum to infinity. That is, we write 



[a (0) = 6 



(E27) 



as a result of [u, T3] = 0. Consequently, the first order 
term for g is, in complete analogy to the discussion lead- 
ing to Eq. 1E12(I . given by 



9 



(i) 



rA. 



(E28) 



mip = mip + TZi 



pi 



(Fl) 



where m; p includes terms in the sum in Eqs. i|D3(l - l|D4|) 

up to a technical cut-off e c while the rest term 1Zi p is the 
sum from s c to infinity computed analytically below. 
At high energies W12 = — W21 ~ 0, while 



Finally, for the second order term g^ 2 \ we have 



[ie n T 3 ,g 



(2)1 



imp 1 y J 



imp ' ; 



(i) 



(E29) 

We solve this equation by using the normalization condi- 
tion, Eq. (|E14|) . Restricting ourselves to isotropic impu- 
rity scattering, we obtain^ 



7(2) 



2e„|eJ 



f 3 iv F ■ VA - A 



i sgn(e„) 



-t 3 



{A-<A) S }),(E30) 



where the inverse scattering time is defined as, 



2ttcN f 



1 + 7T 2 Mlu 2 ' 



(E31) 



For an isotropic (s-wave) superconducting order param- 
eter the last term in Eq. I|E30() vanishes. In this case, the 
second order high-energy contribution from Eq. (|E30|) 
is odd in frequency, and it drops out of the Matsub- 
ara sum. The leading order contribution comes in third 
order— and the high-energy tail of the Matsubara sum is 
oc l/|e„| 3 . This means that the technical cut-off e c can 
be chosen much smaller than in the diffusive limit within 
the Usadel approximation. 

The different high-energy asymptotics within the 
Eilcnbcrger and Usadel approaches are due to the dif- 
fusive approximation employed by Usadel: the impurity 
self-energy, i.e. the inverse scattering time 1/r, is at the 
outset assumed to be the largest energy scale in the prob- 
lem. The high-energy tail is different depending on the 
order in which the limits t — *• and e r — » 00 are taken. 



APPENDIX F: ANALYTIC SUMMATION OF 
THE HIGH-ENERGY TAIL IN THE FOURIER 
SERIES APPROACH 

At high energies e n 3> T c o and J, the matrix W has 
a simple energy dependence that we exploit to sum the 



W11 
W 22 



Ti 



1 + 761 A' 

72 

1 + 7b2A' 



(F2) 
(F3) 



where A 2 = £ n /irT C Q. These relations hold for both the 
trilayer and the pentalayer, which reflects the fact that 
the theory becomes local at high energies (see the effec- 
tive boundary condition ©). The key function of the 
Fourier method then has the form 



3ip 



£s 2 ci + c 2 A 



ds A c 3 + c 4 A + C5A 2 ' 



(F4) 



where 



C2 
C3 
C4 
C5 



7i + (-1)' +P 72 +7i72 [1- 
71762 + (-l) i+p 7276i, 

1 + 71+72 + 7172, 

71762 + 72761 + 761 + 762, 
761762- 



(-1) 



l+p] 



(F5) 
(F6) 
(F7) 
(F8) 
(F9) 



For each element of the matrix m\ v we can perform the 
high-energy Matsubara sum by integration. We get 



5i P -\n[ 1 



c\Jx + c 2 x 



2 1 
—lip, 

s 



TT 2 d 



(F10) 



dx. 



^ ( x + |r) (x + (c 3 + c 4% /i + c 5 x) 



where we used the short hand notation d s = d s /w^s- 
Note that Eq. (|F 10(1 is independent of the temperature 
T and only depends on the parameters in Eqs. l|F5ll -llFC 
on 4 and on the cut-off s r . 
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